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replacing = a + (3 2n+i 'i '^'k~^'~^ (^'('^ ~ 2ra+i )' m=p=2 in the above 
eq. we obtain: 



j=k—n j=k—n 
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Using mean field approximation: {rjkrjj) = {'r]k){'nj) = PkPj = P^', and 
E-;L„^', = 2nPweget: 

j=k—n j=k—n j=k—n 

and given that Pj = = P we obtain the polynomial of 3rd order: 



^ = a - P(a + a' + + ^P' + 4n^(^)P^ 
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Let < a{t) < 1 be the fraction of marked sites. In the mean field descrip- 
tion we can formally define: 



where < oi < 02 < 03... < a2fe+i < 1. 

We choose an odd number of stationary points since we want /(O) > 0, 
/(I) < 0. The odd zeros of f(a), ai,as...a2k+i will be linearly stable fixed 
points while the even number roots will be unstable fixed points. 

If we consider now the effect of DNA replication when the fraction of 
methylated sites is halved then the new "fixed points" corresponding to the 
stable fix points aj will have a fraction of methylated sites right after mitosis 
a* with ao = < a| < (l/2)ai, ....02j < a^j_^i < (l/2)a2j+i, .... 

Let Ti be the period to DNA replication in which the fraction of marked 
sites will increase from a* to 2a* during one cycle then integrating (9), we 
get 
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= U{ai - a{t)) = f{a) 
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where the Bj can be computed in terms of the Oj. 
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In preparing the system we generate a ID lattice of L sites (O's and I's) 
with periodic bounday conditions. For our simulation we used L=1000. 
Each site is a CG nucleotides that has the potential of becoming methy- 
lated or dimethylated. We defined two probability distributions: Pi is the 
probability of site i to be methylated, Pq is the probability of site i to be 
dimethylated. Each probability distribution is constructed based on the 
transition rates of being methylated or dimethylated described in the pa- 
per. The number of neighbors n around each site i that we randomly pick is 
kept fixed, and we sum the number of methylated sites over this neighbors. 
In the simulation we used n=50. In one iteration we've done as following: 
we start with an initial random configuration of sites being methylated and 
dimethylated. We pick up a site i at random, if Yl]tli-n^ij) — -^i '^^^P 
the site methylated, else revert to dimethylated; and if Z]}="_„ < Po 
keep the site dimethylated, else revert to being methylated. I'm doing this 
L times, each time based on previous configuration. I calculate then the 
number of methylated sites and normalize to the length of the lattice. This 
will give me the density of methylated sites at time t. In my simulation 
one unit of time is echivalent with one monte carlo step. To simulate the 
DNA replication process, we are introducing a periodic fluctuation that has 
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as effect the halving of the density of metylated sites periodically at time T. 
In the simulation T is echivalent with 30 monte carlo steps. 
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Abstract 

In aiming to explain the estabhshment, maintenance and stabihty 
of methylation pattern in gene body of Arabidopsis we propose here 
a theoretical framework for understanding how the methylated and 
unmethylated states of cytosine residues are maintained and trans- 
mitted during DNA replication. Routed in statistical mechanics, the 
framework built herein is used to explore minimal models of epigenetic 
inheritance and identify the necessary conditions for stability of methy- 
lated/unmethylated states of cytosine over rounds of DNA replication. 
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The models are flexible enough to allow adding new biological concepts 
and information. 

Keywords: lattice models, mean field analytical approach, Monte Carlo 
simulations. 

1 Introduction 

DNA methylation is a dynamic epigenetic process that refers to the en- 
zymatic transfer of a methyl group to the specific nucleotides within the 
DNA sequence. In eukaryotes, this modification or marks affects almost ex- 
clusively cytosines [1]. DNA methylation is readily detected in plants and 
mammals, where it is critical for normal development and genome stability. 
Interestingly, plants seem more prone to the inheritance of DNA methylation 
defects than mammals [2,3]. 

Due to a near-complete genome sequence annotated to very high stan- 
dards, a comprehensive set of genomics tools and powerful genetics, the 
fiowering plant Arabidopsis has rapidly become a prime model for the study 
of DNA methylation and its inheritance patterns in higher eukaryotes. In 
this work we will refer at Arabidopsis as a model organism although the 
framework in which we work it can be applied to DNA methylation in other 
organisms as well. 
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In the Arabidopsis, methylation of cytosine has been detected on gene 
body, gene promoters and repeat elements (transposable elements). If the 
role of methylation in the context of repeat elements is considered to be 
of defense against invasive DNA and on gene promoters of silencing the 
gene, the role of methylation on gene body is not yet clear. On gene body 
methylation of cytosine is restricted to CG sites, [4,5] by difference with 
the repeat elements or gene promoters case where methylation is sequence 
dependent and can be found also on CHG and CHH sites (where H can be 
any of the four nucleotides: A,T,C,G) [6,7] The relative prevalence of DNA 
methylation in each sequence context throughout the genome was assessed, 
revealing that 55% were in CG context, while 23% and 22% were in the 
CHG and CHH contexts, respectively [4,5,8,9]. 

Gene body methylation occurs on about a third of all genes, and these 
genes tend to be highly and ubiquitously expressed in different Arabidopsis 
tissues [10,11] 

One of the defining properties of epigenetic phenomena is its stability, 
the ability of the cell to maintain its epigenetic state stable through many 
cell divisions. The density of the marks, (methylated states) responsible 
for the epigenetic effects, is changed during DNA replication by introducing 
newly synthesized DNA indicating that these heritable states must be robust 
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against significant perturbations in there concentration. In the same time, 
mechanisms of DNA methylation involve enzymes that can act on more than 
one nucleotide in its neighborhood. This non-locality of action opens the 
possibility of interesting collective aspects that have a role in maintaining 
the stability of epigenetic states. 

We approach the problem by using methods that traditionally are used in 
statistical mechanics and dynamical systems. In my previous work on epige- 
netic processes [12], I was mostly concerned with understanding the stability 
of histone modifications, rather than DNA methylation. Given that DNA 
methylation is another important epigenetic mark critical in development 
and genome stability, in this work we wish to explore the stability of DNA 
methylation pattern across multiple generations and focus on understanding 
DNA methylation in a first approximation in the context of gene body. 

The aim is to explore the properties of a minimal model of epigenetic 
silencing in order to identify the necessary conditions for stability of cytosin 
states that correspond to distinct epigenetic phenotypes i.e. methylated/ 
un-methylated states. The model is based on the current understanding of 
DNA methylation in gene body Arabidopsis, with particular emphasis on 
the interplay between the mechanisms that enable the establishment and 
maintenance of this modifications [13]. 
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In section 2 we will present the general framework of the model. In 
order to make the present discussion self-contained we will introduce some 
well known aspects of the methods commonly used in statistical physics (see 
also [12]). In section 3 we will apply the framework and methods presented in 
previous section to the context of gene body methylation in Arabidopsis and 
present the results. In section 4 we will discuss some aspects of the model 
and in section 5 we will draw the conclusions and present some posible future 
directions of the model. 

2 Methods and general framework of the model 

We consider a ID lattice of size L whose sites correspond to nucleotides/cytosines 
ordered along the length of the DNA. The nucleotide corresponding to site i, 
can be in several states, corresponding to particular situation that we are in- 
terested in. These states are labeled by s = 1, iV. The rates of transition 
at site i from state s' to state s, namely, Riss>{si, . . . , s^-i, s', Sj+i, . . . , s^], 
depends not only on the local state but also on the states of all the neighbors 
within a range I. In practice, this dependency arises because particular mod- 
ifications of a site leads to recruitment of particular enzymes that could affect 
modification rates of the neighboring nucleotides. The master equation de- 
scribing the time evolution of the probability distribution P[si, . . . ,SL]t\ is 
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given by: 

j^P[su...,SL;t] = 

L 

= {^Sis'[si, ■ ■ . ,Si-i,s',Si+i, . . .,Sl]P[si, . . . ,Si-i,s',Si+i, ...,SL;t\ 

i=l s' 

-Ris'siisi, Si-i,Si, Si+i, sl]P[si, Si-i,Si, Si+i, ...,sl; t]) (1) 

for times between DNA replication. At the point of DNA duplication, a 
novo strand is formed. This will have as consequences that the fraction 
of methylated sites right after DNA replication will be diluted. Taking in 
consideration that right after DNA replication we have a hemimethylated 
DNA, we represent that in the evolution of probability distribution. Actually 
for the gene body case, density of methylated sites will be halved, i.e will be 
half the fraction of methylated sites that were before DNA replication. In 
this process we assume that DNA duplication happens instantaneously (in 
reality, fast compared to the time between two duplication events). Fig. 1 
provides a schematic representation of the model and its dynamics. 

To solve the master equation (1) analytically for the long time behavior of 
P[si, . . . , sl; t] is generally an impossible task. One, therefore, has to resort 
to some sort of approximation. One such approximation often used success- 
fully in statistical mechanics is the "mean field" approximation [14]. In this 
approach one approximates P[si, . . . , SL',t] by a f actor ized form niPil^i! 
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Figure 1: Methylation on CG sequence is symetrically on both strands; right after 
DNA replication at T+ we have a hemimethylated DNA and density of methylated 
sites is halved; L - length of the lattice; 1 - range of cooperative behavior; T - 
periodic time to DNA replication 

Using this approximation one derives that the evolution equation for Pi[si; t] 
is going to be 

^Pi[si;t] = ^{Ris,s'Pi[s']t\ - Ris's,Pi[si]t]) (2) 

s' 

where the definition of the average rates Riss' is 

Riss' = ^ Riss'[si, . . . ,SL]pi[si;t] . . .pi-i[Si^i]t]pi+i[si+i,t] . . .pL[sL;t]. 

Sl,...,Si_l,Si + l,...,S£, 

(3) 

Notice that these averaged rates Riss' are polynomials in making 
eq.(2) a nonlinear equation. 
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In the mean field analysis of all the models discussed in this work, we 
will ignore the spatial variation of 'marks' and replace them by average 
concentrations corresponding to an entire region of DNA, namely Pi[si; t] = 
p[si;t]. We thereby focus on regions of DNA with one epigenetic fate and 
be concerned with 'uniform' states. The equations for the variables p[s; t] 
are: 

j^p[s; t] = Y,i^ss'P[s'; t] - Rs'sP[s; t]) (4) 

s' 

where Rss' = Riss', is given by Eq. [2l These are independent of i because 
the rules of transitions are translation invariant and we ignore boundary 
effects. On incorporating recruitment and cooperative behavior multiple 
neighboring sites of a site influence the probability of the state at that site, 
therefore, the transition rates are dependent on what happens on neigh- 
boring sites. We suppose that the rates Ris^s'isi, ■ ■ ■ , s', Sj+i, . . . ,sl] 
depend only on the fraction of sites in a given state in the neighborhood 
of i within separation I, where 1 « I (we could still have / << L to be 
physically meaningful). We can group then L sites into L/l clusters of / sites 
each, i.e. coarse-graining the system. We redefine the probabilities pi[si,t] 
of state Si at site z G [1; by the averaged probability Pj[s, t] of state S at 
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any cluster j G [l,L/l], where formally 

1 

i=jl-l+l 

Moreover we can assume that the averaged probabilities are approximately 
site independent. The new states S are not binary corresponding to the 
presence or absence of marks but a discrete spectrum of states that can be 
approximated by the concentration of marks in a cluster. This mean-field 
equivalence of the local probability of a binary state at a site to the prob- 
ability density (or normalized concentration) of states in a 'coarse-grained 
cluster' is going to be exploited in the rest of the work implicitly in writing 
down mean-field differential equations for the dynamics of the system. We 
will not introduce in the rest of the work the formal redefinitions of proba- 
bilities done above. The mean field approximation, turns out, a posteriori, 
to be justified and quite effective in many cases [14,15]. This method as 
shown is based on averaged quantities that coarse -grain the system and by 
neglecting the spontaneous fluctuations in the concentration of the states, 
predicts long range order. 

We will study, analytically and computationally the stochastic model 
of epigenetic inheritance formulated above for a particular choice of states 
and rules of state transitions proper for describing DNA methylation in 
gene body Arabidopsis. In the next section the discussion will be on a 
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concrete case of DNA methylation in gene body Arabidopsis. Here we will 
show that some restrictions in the dynamics (transition rates) and some 
addtional constraints are required for the recovery of the epigenetic marks 
to take place. 

3 Results: Modeling gene body methylation in 
Arabidopsis 

Abiding by our goal of identifying minimal models of epigenetic DNA methy- 
lation, we develop in this section a two-state model for studying stable epige- 
netic marks and understanding gene body methylation. In Arabidopsis gene 
body, methylation is restricted almost exclusively to CG sites and seams to 
be associated with expression rather then silencing [4,5]. 

Following the understanding of Colot group et. al.[16] for DNA methy- 
lation in Arabidopsis we are considering that the process of methylation 
takes places in two critical steps. First step concerns establishment of DNA 
methylation pattern and its associated mechanisms; while the second con- 
cerns maintenance of this modification within and between generations. 

Based on the general framework presented in previous section we consider 
the string of nucleotides as a 1 dimensional finite lattice that approximate 
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the DNA. As experimentally has been established that methylation status 
is influenced by nearby cytosine [17] we have to take into account this fact 
in generating the model of dynamic evolution of the system. Thus, in our 
model, the rates of transition at site k depends of the states of all the 
neighbors cytosine within a range l=n, property called cooperativity and 
of an inherent constant defined by the enzymes involved in the establishing 
mechanism of DNA methylation. The inherent constant defines the property 
of a cytosine to become methylated or de-methylated independently of its 
neighbors. This can be considered a de novo methylation where a cytosine is 
methylated/de- methylated with a constant rate due to enzymatic machinery. 
By contrast the cooperative term defines the dependency of the state of the 
cytosine at site k of the methylated status of its n neighbors. Such a term 
describes the local modulation of the enzymatic machinery by raising or 
lowering the local concentration of enzymes at a given place in the lattice. 

In this sense the dynamics that determine the establishment of DNA 
methylation, the transition rate at site k is defined as bellow by the two 
components, see fig. 2: the component that defines the inherent property 
of cytosine to become methylated or un-methylated independently of its 
neighbors that we call INHERENT rate and COOPERATIVITY component 
which is determined and therefore depends by the state of the n neighbors 
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that surround the cytosine, 

cooperativity 
inherent ^ ^ 

Mathematically all of above are written as following: = a +/3 ; 

2n + 1 



r_k=inherGnt rate m+coop rate m 




r'_k=inherent rate_u+coop rate_u 



Figure 2: transition rates at site k; m-methylated state; u-unmethylated state; 
is the rate of unmethylated cytosine at site k to become methylated; rj. the rate of 
methylated cytosine at site k to become unmethylated 

rifc = Yl^t.k-n ^i(^)' ^ is inherent rate constant of cytosines to be methy- 
lated; /3 is a proportionality constant; rj is methylated state and (1 — rj) 
-de-methylated state. 

r'^ = a' + I3'{1 - 2^); nk = I]jifc-n^i(*); depends of n, number of 
neighbors; where a' is the inherent constant rate of de-methylation ; /3' is a 
proportionality constant of de-methylation. 

Given the dynamics described above we write the equation for time evo- 
lution of the density of methylated sites for times between DNA replication 



12 



known as master equation; with (rjk) density of methylated sites: 



%^ = ((l-%)r,)-(%r^); 
^ = ((1 - Vk){a + P^^)) - iVkla ^^)]); 

1 / \ n h-\-Ti Q k-\--Ti Qf k-\-n 

j=k—l j=k—l j=k—l 

(6) 

To solve analytically equation (6) we use mean field approximation men- 
tioned in the previous section, where = = PkPj,k / j and 
also Pj = Pk = P; to obtain a simpler equation that describes the dynamics 
of density of methylated sites at times between DNA replication: 

dt 2n + l ^ ^ ^ 2n + l ^ j A; ^ ^/^ ; 2n + 1 ^ ^ 

j=fc— 1 j=fc— 1 j=k—l 

or 

^ = a + 2nP - aP ^2nP2 - (a' + /3')P + ^^2nP2; 

2n + 1 2n + 1 ^ ^ ^ 2n + 1 

and by grouping the terms: 

V ' V ' 

In steady state this model is a quadratic equation therefore at it's best can 
have just one stable state even in the absence of fluctuations induced by 
perturbations due to DNA replication. 

13 



Epigenetic DNA methylation implies alternative states that are stable 
over time and are inherited through cell division. Any model that tries to 
explain the methylation process should be able to obtain a coexistence of 
stable states. Actually the understanding of epigenetic processes, in terms of 
multiple steady states, has been suggested already long ago by Waddington 
and most clearly by Delbrucks [19,20] Multistationarity is the property of 
systems whose structure is such that they can display two or more distinct 
steady states under identical conditions. In our model the requirement of 
mulitistationarity (or at least of a bistable state) can happen if we allow 
the de-methylated and methylated sites to recruit enzymes cooperatively 
in a non-linear manner to de-methylated and methylated neighboring sites 
respectively. 

This will affect the transition rates for methylation/ de-methylation 
to include a degree of non-linear cooperative methylation, respectively de- 
methylation: = a + P^^; r'^ = a' + /3'(1 - nu = E^tk-nVjit) 

where m is the degree of non-linear cooperative methylation and p is the 
degree of non-linear cooperative de-methylation. 

In this new context, the quadratic equation (7) is now changed to a 
polynomial of higher order and the model is modified to enable the presence 
of multiple dynamical attractors. Even for the simplest case of non-linear 
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cooperative behavior when m=p=2 we get in steady state a polynomial of 
order 3 instead of the quadratic equation (7) (see Appendix A), obtaining 
a model that can have a bi-stable state and as such the requirement for 
epigentic DNA methylation and memory. 

In fact if we are calling f(a) the right side of such polynomial equation 
for m=p=2, f(a) will have three zeros, ai < a2 < in the interval (0,1). 
The scenario relevant to us is when ai and 03 are stable and are separated 
by 02 unstable (03 corresponds to high concentration of marks and oi to low 
concentration of marks) see fig. 3 
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Figure 3: polynomial equation of order 3 when m=p=2 

Any initial states with a(0) < 02 will eventually be attracted to ai while 
any initial state with a(0) > 02 will be attracted to a-^. Suppose now that 
DNA replication takes place periodically at time T. Right after replication 
the density of methylated sites is halved (see fig.l). It is clear then that if 
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«2 > 2*^3' tlisre will be only one stable fixed point, which will be close to ai 
if T is large enough (time T to replication considerably larger then the time 
scale of methylation rates). To see this we simply note that starting from 
as the cell after DNA replication will have a value which is less than 02 and 
so will enter the basis of atraction of the stable fixed point ai. 

However for 02 < ^ and T fulfilling the same conditioned stated earlier, 
then there will be two stable fixed point, one near oi and one near 03 even 
with perturbations induced by DNA replication. These puts restrictions on 
the parameters entering f(a) and T required for the existence of multiple 
stable points and thus of epigenetic memory and regulation; the density of 
methylated sites right after DNA replication (after halving) has to have a 
higher value then 02 the unstable value. We shall not go into this here. It 
is also possible to give fairly explicit expressions for T in terms of f (a) , (see 
Appendix B) 

Going beyound mean field, using Monte Carlo simulations we would like 
to see if the requirement of stability is mainteined in the above framework 
when the perturbations induced by DNA replications are taken in consid- 
eration. In this sense fig 4. shows how methylation pattern is established 
and remain stable, reaching a stationary state from following the dynamics 
introduced so far. The DNA replication is not yet involved in the process. 
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establishement of methylated pattern: evolution of density of methylated states 
0.65 1 , , 1 , 1 




time 



Figure 4: Establishment of metliylated pattern in the absence of DNA replication- 
evolution of the density of methylated states; lunit time « 1 mcs. 

When we introduced in the system the perturbation due to DNA reph- 
cation by simply halving periodically the density of methylated sites we see 
in fig. 5 that the methylated pattern, the density of methylated sites is 
recovered, provided that the time to replication is longer then the recovery 
time, showing as such, explicitly that the dynamics described suffix for the 
stability of the methylated pattern. 

After multiple replication cycles one expects methylation to not be ex- 
actly as it was originally because of accidental loss or gain, see fig. 5. For 
Monte Carlo simulation details and how the system was prepared see Ap- 
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maintenance of of metfiylated pattern: periodic recovery of methylated states 
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Figure 5: maintenance of methylated pattern in the presence of DNA replication- 
periodic recovery of methylated states; lunit time ~ 1 mcs. 

pendix C. All the parameters used in the simulation were chosen based on 
literature research. 



4 Discussions 



The requirement of non-linearity to can have multiple dynamical attractors 
in the system has been previously studied in other epigenetic contexts [21]. 
This studies showed that to obtain multistationarity in the system a pos- 
itive feedback circuit is a necessary condition. We note, however, that at 
a fundamental level the presence of non-linear cooperative recruitment of 
enzymes to methylated/de-methylated states might be due to the existence 
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of such positive feedback circuit. In Arabidopsis more states are introduced 
when methylation takes place on repeat elements over their entire length 
and strongly correlates with transcription inhibition. Methylation in this 
case is on CG, CHG and CHH sites instead of solely CG sites like on gene 
body, and as such is sequence dependent. We know from previous work [12] 
that the presence of intermediate states naturally lead to cooperative effects 
when each of the intermediate states recruit enzymes for further modifica- 
tion. A continuation of this work would be to analyse the effects of sequene 
dependency on the stability of the methylation pattern and try to under- 
stand if non-linearity in the transition rates is still a key ingredient in the 
maintenance of DNA methylation or if solely the presence of more states is 
stabilizing the system. We have phrased the mean-field theory in terms of 
coarse grained quantities like the fraction of sites with a particular mark in 
a cluster. Given that methylation on repeat element takes place over their 
entire length, a natural question would be to try understand how does the 
effective model change if we continue the coarse-graining to larger length 
scales. In other words: one could ask how the model renormalizes under 
iterative blocking transformations [14,15]. In absence of any conservation 
law, there is no obvious reason why this system should not have a finite 
(although long) correlation length in space and, similarly, a finite correla- 
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tion time. The system would not have genuinely multiple phases. All these 
effects, which are missed by mean-field theory, would, in principle, show up 
in renormalization group. 

5 Conclusions 

The stability of epigenetic DNA methylation is a rich subject in biology. 
The exact function(s), of much of the DNA methylation found outside of 
repeat elements remain unclear [6,7]. Using an approach routed in statis- 
tical physics we proposed a theoretical framework for understanding how 
the methylated and non-methylated states of cytosine are maintained and 
transmitted when perturbations (as DNA replication) are involved. 

The model explains the establishment of DNA methylation pattern i.e. 
study the dynamics of the density of the methylated sites. Analyzes the ef- 
fect of the non-linear cooperativity (in the transition rates) on the stability 
of the marks and shows that at least in the gene body case where methy- 
lation is restricted to CG sites in order to have stable DNA methylation 
patterns transmitted over generations non-linear cooperativity is required 
in the maintenance process. 

Many features presented here in the context of Arbidopsis can be also 
extendet to DNA methylation in other organisms as well. The model ex- 
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tends the view that multi-stationarity in gene body DNA methylation pat- 
tern arises by aUowing the de-methylated and methylated sites to recruit 
enzymes cooperatively in a non-linear manner to de-methylated and methy- 
lated neighboring sites respectively. We finish by concluding that overall 
the present work lays the grounds for understanding DNA methylation on 
repeat elements in Arabidopsis, and as such extending the mathematical 
framework to compleate the modelling and understanding of DNA methy- 
lation in Arabidopsis. 
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